Go to the end of this file for the parameter estimates of the best fitted model on phenotype width
Trying to fit non-linear model with BRMS(Bayesian Regression Models using Stan)
widths <- read_csv("widths.csv")
widths
sample.list <- unique(widths$ID)
lapply(1:6, function(i) {
ggplot(widths) +
geom_line(aes(x=days, y=width, group=ID)) +
facet_wrap_paginate(~ID, ncol = 6, nrow = 5, page =i)
}
)
First set up the formula
weibull.bf1 <- bf(
width ~ Hmax - (Hmax - Hmin) * exp(-(k*days)^delta), #Weibull model
Hmax + Hmin + k + delta ~ 1, #general model, paramters do not vary for individuals
nl=TRUE)
Priors. Hmin and Hmax use median at start and end dates.
stat.data <- function(data) {
data %>%
group_by(days) %>%
summarize(meadian=median(width),
max=max(width),
min=min(width),
sd=sd(width))
}
stat.data(widths)
summary(fit1, waic=TRUE, R2=TRUE)
So with these priors it seems that the model can flip between having high and low Hmin and Hmax. This relates to the delta parameter flipping signs. So one possibility is to keep delta positive. Or more tightly constrain the priors on Hmax and Hmin.
summary(fit2, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-(k * days)^delta)
Hmax ~ 1
Hmin ~ 1
k ~ 1
delta ~ 1
Data: widths (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 5042.2; R2 = 0.32
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 59.78 0.71 58.43 61.17 344 1.00
Hmin_Intercept 29.69 5.14 16.08 36.98 50 1.10
k_Intercept 0.01 0.00 0.01 0.01 64 1.08
delta_Intercept 4.30 0.74 2.90 5.81 64 1.08
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 10.78 0.29 10.21 11.34 170 1.05
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit2)
pairs(fit2)
rescale k to get it more reasonable
weibull.bf2 <- bf(
width ~ Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta), #Weibull model
Hmax + Hmin + k + delta ~ 1, #general model, paramters do not vary for individuals
nl=TRUE)
summary(fit3, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta)
Hmax ~ 1
Hmin ~ 1
k ~ 1
delta ~ 1
Data: widths (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 5040.41; R2 = 0.33
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 59.67 0.64 58.45 60.97 2415 1.00
Hmin_Intercept 31.18 5.50 17.27 38.14 692 1.00
k_Intercept 0.81 0.05 0.75 0.93 727 1.01
delta_Intercept 4.86 0.95 3.10 6.73 988 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 10.74 0.30 10.17 11.35 2360 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit3)
pairs(fit3)
What does the fit look like?
newdata <- data.frame(days=seq(min(widths$days), max(widths$days),1))
fit3.fitted <- cbind(newdata, fitted(fit3, newdata)) %>% as.tibble() %>%
rename(width=Estimate, lower.ci='2.5%ile', upper.ci='97.5%ile')
plot
pl <- ggplot(aes(x=days, y=width),data=NULL)
pl <- pl + geom_line(aes(group=ID), alpha=.2, data=widths)
pl + geom_line(color="skyblue", lwd=1.5, data=fit3.fitted)
Now try adding random effects for model parameters:
What parameters do we think might be interesting to allow to vary? Probably not Hmin. Try making a series of plots to see how varying delta or k affects things:
weibull.fn <- function(Hmax,Hmin,k,delta,days) {
Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta) #Weibull model
}
for(delta1 in seq(3,5,.25)) {
tmp.widths <- weibull.fn(Hmax=62,
Hmin=20,
k=.7,
delta=delta1,
days=newdata$days)
abc <- data.frame(newdata$days, tmp.widths)
p <- ggplot(data=abc, aes(newdata$days, tmp.widths)) +
geom_line() + ylim(0,100) + ggtitle("delta=",delta1)
print(p)
}
for(k1 in seq(.5,1.5,.1)) {
tmp.widths <- weibull.fn(Hmax=60,
Hmin=20,
k=k1,
delta=4,
days=newdata$days)
abc <- data.frame(newdata$days, tmp.widths)
p <- ggplot(data=abc, aes(newdata$days, tmp.widths)) +
geom_line() + ylim(0,80) + ggtitle("k=",k1)
print(p)
}
First try with only fixing Hmin
weibull.bf3 <- bf(
width ~ Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta), #Weibull model
Hmax + k + delta ~ (1|ID), # vary for individuals
Hmin ~ 1, # do not vary per individual
nl=TRUE)
summary(fit5, waic=TRUE, R2=TRUE)
There were 12 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gaussian(identity)
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta)
Hmax ~ (1 | ID)
k ~ (1 | ID)
delta ~ (1 | ID)
Hmin ~ 1
Data: widths (Number of observations: 664)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
ICs: LOO = NA; WAIC = 4943.66; R2 = 0.54
Group-Level Effects:
~ID (Number of levels: 166)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept) 3.94 1.07 1.34 5.81 526 1.01
sd(k_Intercept) 0.14 0.02 0.10 0.17 1611 1.00
sd(delta_Intercept) 0.64 0.59 0.02 2.22 4596 1.00
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 59.81 0.68 58.44 61.12 2333 1.00
k_Intercept 0.83 0.02 0.79 0.87 2558 1.00
delta_Intercept 7.54 0.51 6.52 8.52 10000 1.00
Hmin_Intercept 33.01 1.58 29.81 35.95 4256 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 8.84 0.33 8.23 9.51 2920 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
get fitted values
b <- widths %>% ungroup()
fit5.fitted <- cbind(b, fitted(fit5)) %>% as.tibble()
fit5.fitted
plot
plot.fitted <- function(fitted.data) {
pl <- ggplot(fitted.data, aes(x=days)) +
geom_line(aes(y=width),color="blue") +
geom_line(aes(y=Estimate),color="red") +
facet_wrap_paginate(~ID, ncol = 6, nrow = 5)
pages <- n_pages(pl)
lapply(seq_len(pages),
function(i) {
ggplot(fitted.data, aes(x=days)) +
geom_line(aes(y=width), color="blue") +
geom_line(aes(y=Estimate), color="red") +
facet_wrap_paginate(~ID, ncol = 6, nrow = 5, page=i)
}
)
}
plot.fitted(fit5.fitted)
plot fitted vs actual:
plot.fitted.actual <- function(fn) {
r_squared <- summary(lm(Estimate ~ width, data=fn))$adj.r.squared
r_squared <- round(r_squared, digits=4)
fn %>%
mutate(days=as.factor(days)) %>%
ggplot(aes(x=width, y=Estimate, shape=days, color=days)) +
geom_point() +
geom_abline(intercept = 0, slope=1) +
ggtitle(paste0("R2 = ",r_squared))
}
plot.fitted.actual(fit5.fitted)
total_sum_of_squared_residuals <- function(fn) {anova(lm(Estimate ~ width, data=fn))[2,2]}
cv
TSSR.fit <- total_sum_of_squared_residuals(fit5.fitted)
waic.fit <- round(waic(fit5)$waic, digits=2)
kfoldic.fit <- round(kfold(fit5)$kfoldic, digits=2)
widths.mean <- read_csv("widths.mean.csv")
widths.mean
lapply(1:6, function(i) {
ggplot(aes(x=days, group=ID), data=widths.mean) +
geom_line(aes(y=width), color="red") +
geom_line(aes(y=width_raw), color="black") +
facet_wrap_paginate(~ID, ncol = 6, nrow = 5, page =i)
}
)
stat.data(widths.mean)
stat.data(widths)
summary(fit3.mean, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta)
Hmax ~ 1
Hmin ~ 1
k ~ 1
delta ~ 1
Data: widths.mean (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 4913.98; R2 = 0.48
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 65.46 0.84 63.98 67.26 2172 1.00
Hmin_Intercept 31.96 3.81 22.55 37.64 1081 1.00
k_Intercept 0.77 0.03 0.73 0.83 1355 1.00
delta_Intercept 4.60 0.77 3.20 6.16 1310 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 9.76 0.27 9.26 10.31 2722 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit3.mean)
pairs(fit3.mean)
What does the fit look like?
newdata <- data.frame(days=seq(min(widths.mean$days), max(widths.mean$days),1))
fit3.mean.fitted <- cbind(newdata, fitted(fit3.mean, newdata)) %>%
as.tibble() %>%
rename(width=Estimate, lower.ci='2.5%ile', upper.ci='97.5%ile')
plot
pl <- ggplot(aes(x=days, y=width),data=NULL)
pl <- pl + geom_line(aes(group=ID), alpha=.1, data=widths.mean)
pl + geom_line(color="skyblue", lwd=1.5, data=fit3.mean.fitted) +
ggtitle("width.mean vs. days")
Now try adding random effects for model parameters:
What parameters do we think might be interesting to allow to vary? Probably not Hmin. Try making a series of plots to see how varying delta or k affects things:
weibull.fn <- function(Hmax,Hmin,k,delta,days) {
Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta) #Weibull model
}
for(delta1 in seq(3,8,.5)) {
tmp.widths <- weibull.fn(Hmax=65,
Hmin=40,
k=1,
delta=delta1,
days=newdata$days)
abc <- data.frame(newdata$days, tmp.widths)
p <- ggplot(data=abc, aes(newdata$days, tmp.widths)) +
geom_line() + ylim(0,80) + ggtitle("delta=",delta1)
print(p)
}
for(k1 in seq(0,5,.5)) {
tmp.width <- weibull.fn(Hmax=65,
Hmin=40,
k=k1,
delta=5,
days=newdata$days)
abc <- data.frame(newdata$days, tmp.widths)
p <- ggplot(data=abc, aes(newdata$days, tmp.widths)) +
geom_line() + ylim(0,80) + ggtitle("k=",k1)
print(p)
}
try with only fixing Hmin
weibull.bf3 <- bf(
width ~ Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta), #Weibull model
Hmax + k + delta ~ (1|ID), # vary for individuals
Hmin ~ 1, # do not vary per individual
nl=TRUE)
summary(fit5.mean, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta)
Hmax ~ (1 | ID)
k ~ (1 | ID)
delta ~ (1 | ID)
Hmin ~ 1
Data: widths.mean (Number of observations: 664)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
ICs: LOO = NA; WAIC = 3994.06; R2 = 0.92
Group-Level Effects:
~ID (Number of levels: 166)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept) 7.80 0.61 6.64 9.04 1874 1.00
sd(k_Intercept) 0.13 0.01 0.11 0.15 2589 1.00
sd(delta_Intercept) 2.49 0.36 1.85 3.25 1669 1.00
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 65.89 0.76 64.38 67.37 2669 1.00
k_Intercept 0.85 0.01 0.82 0.88 2056 1.00
delta_Intercept 5.43 0.42 4.63 6.28 4511 1.00
Hmin_Intercept 23.27 1.41 20.34 25.83 5660 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 3.78 0.22 3.38 4.23 993 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
get fitted values
b <- widths.mean %>% ungroup()
fit5.mean.fitted <- cbind(b, fitted(fit5.mean)) %>% as.tibble()
fit5.mean.fitted
plot
plot.fitted(fit5.mean.fitted)
plot fitted vs actual(modified):
plot.fitted.actual(fit5.mean.fitted)
plot fitted vs actual(raw):
plot.fitted.actual.raw <- function(fn) {
r_squared <- summary(lm(Estimate ~ width_raw, data=fn))$adj.r.squared
r_squared <- round(r_squared, digits=4)
fn %>%
mutate(days=as.factor(days)) %>%
ggplot(aes(x=width_raw, y=Estimate, shape=days, color=days)) +
geom_point() +
geom_abline(intercept = 0, slope=1) +
ggtitle(paste0("R2 = ",r_squared))
}
plot.fitted.actual.raw(fit5.mean.fitted)
TSSR on raw data set
total_sum_of_squared_residuals_raw <- function(fn) {anova(lm(Estimate ~ width_raw, data=fn))[2,2]}
TSSR.fit.mean.raw <- total_sum_of_squared_residuals_raw(fit5.mean.fitted)
cv
TSSR.fit.mean <- total_sum_of_squared_residuals(fit5.mean.fitted)
waic.fit.mean <- round(waic(fit5.mean)$waic, digits=2)
kfoldic.fit.mean <- round(kfold(fit5.mean)$kfoldic, digits=2)
widths.pmm <- read_csv("widths.pmm.csv")
widths.pmm
lapply(1:6, function(i) {
ggplot(aes(x=days, group=ID), data=widths.pmm) +
geom_line(aes(y=width), color="red") + #red: modified data
geom_line(aes(y=width_raw), color="black") +
facet_wrap_paginate(~ID, ncol = 6, nrow = 5, page =i)
}
)
Priors. Hmin and Hmax use median at start and end dates.
stat.data(widths)
stat.data(widths.mean)
stat.data(widths.pmm)
the stats of widths.mean and widths.pmm are very similar. so start the prior of width.mean.
weibull.bf3 <- bf(
width ~ Hmax - (Hmax - Hmin) * exp(-((k/100)*days)^delta), #Weibull model
Hmax + k + delta ~ (1|ID), # vary for individuals
Hmin ~ 1, # do not vary per individual
nl=TRUE)
summary(fit5.pmm, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: width ~ Hmax - (Hmax - Hmin) * exp(-((k/100) * days)^delta)
Hmax ~ (1 | ID)
k ~ (1 | ID)
delta ~ (1 | ID)
Hmin ~ 1
Data: widths.pmm (Number of observations: 664)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
ICs: LOO = NA; WAIC = 4269.48; R2 = 0.88
Group-Level Effects:
~ID (Number of levels: 166)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept) 7.31 0.64 6.11 8.61 1578 1.00
sd(k_Intercept) 0.13 0.01 0.11 0.15 2101 1.00
sd(delta_Intercept) 2.43 0.38 1.74 3.21 1968 1.00
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 65.43 0.74 63.95 66.87 2204 1.00
k_Intercept 0.84 0.02 0.81 0.87 1922 1.00
delta_Intercept 5.84 0.46 4.97 6.78 4206 1.00
Hmin_Intercept 25.49 1.53 22.27 28.33 3241 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 4.76 0.25 4.29 5.27 1593 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
fitted values
a <- widths.pmm %>% ungroup()
fit5.pmm.fitted <- cbind(a, fitted(fit5.pmm)) %>% as.tibble()
plot fitted vs actual(modified)
plot.fitted.actual(fit5.pmm.fitted)
plot fitted vs actual(raw)
plot.fitted.actual.raw(fit5.pmm.fitted)
CV
TSSR.fit.pmm.raw <- total_sum_of_squared_residuals_raw(fit5.pmm.fitted)
TSSR.fit.pmm <- total_sum_of_squared_residuals(fit5.pmm.fitted)
waic.fit.pmm <- round(waic(fit5.pmm)$waic, digits=2)
#### which fitted model is the better?
create an empty list for models
cv.list[2,1] <- TSSR.fit.mean.raw
cv.list[3,1] <- TSSR.fit.pmm.raw
cv.list[1,2] <- TSSR.fit
cv.list[2,2] <- TSSR.fit.mean
cv.list[3,2] <- TSSR.fit.pmm
WAIC(Widely Applicable Information Criterion) an extension of the Akaike Information Criterion (AIC) that is more fully Bayesian.
K-fold CV Data are randomly partitioned into K subsets of equal size. Then the model is refit 10 times(default), each time leaving out one of the 10 subsets.
cv.list[1,3] <- waic.fit
cv.list[2,3] <- waic.fit.mean
cv.list[3,3] <- waic.fit.pmm
cv.list[1,4] <- kfoldic.fit
cv.list[2,4] <- kfoldic.fit.mean
cv.list[3,4] <- kfoldic.fit.pmm
cv.list
#### parameter estimates of best fitted model on phenotype width : take fit5.mean
read
1) population (fixed) & group level (random) effects
2) parameters: Hmax, k & delta of each genotype
width.best.fitted <- fit5.mean$fit
dimnames(width.best.fitted)
$iterations
NULL
$chains
[1] "chain:1" "chain:2" "chain:3" "chain:4"
$parameters
[1] "b_Hmax_Intercept" "b_k_Intercept" "b_delta_Intercept"
[4] "b_Hmin_Intercept" "sd_ID__Hmax_Intercept" "sd_ID__k_Intercept"
[7] "sd_ID__delta_Intercept" "sigma" "r_ID__Hmax[ID_1,Intercept]"
[10] "r_ID__Hmax[ID_10,Intercept]" "r_ID__Hmax[ID_100,Intercept]" "r_ID__Hmax[ID_101,Intercept]"
[13] "r_ID__Hmax[ID_102,Intercept]" "r_ID__Hmax[ID_103,Intercept]" "r_ID__Hmax[ID_104,Intercept]"
[16] "r_ID__Hmax[ID_105,Intercept]" "r_ID__Hmax[ID_106,Intercept]" "r_ID__Hmax[ID_107,Intercept]"
[19] "r_ID__Hmax[ID_109,Intercept]" "r_ID__Hmax[ID_111,Intercept]" "r_ID__Hmax[ID_112,Intercept]"
[22] "r_ID__Hmax[ID_113,Intercept]" "r_ID__Hmax[ID_114,Intercept]" "r_ID__Hmax[ID_116,Intercept]"
[25] "r_ID__Hmax[ID_117,Intercept]" "r_ID__Hmax[ID_118,Intercept]" "r_ID__Hmax[ID_119,Intercept]"
[28] "r_ID__Hmax[ID_12,Intercept]" "r_ID__Hmax[ID_120,Intercept]" "r_ID__Hmax[ID_121,Intercept]"
[31] "r_ID__Hmax[ID_122,Intercept]" "r_ID__Hmax[ID_123,Intercept]" "r_ID__Hmax[ID_124,Intercept]"
[34] "r_ID__Hmax[ID_125,Intercept]" "r_ID__Hmax[ID_126,Intercept]" "r_ID__Hmax[ID_127,Intercept]"
[37] "r_ID__Hmax[ID_128,Intercept]" "r_ID__Hmax[ID_129,Intercept]" "r_ID__Hmax[ID_13,Intercept]"
[40] "r_ID__Hmax[ID_130,Intercept]" "r_ID__Hmax[ID_132,Intercept]" "r_ID__Hmax[ID_133,Intercept]"
[43] "r_ID__Hmax[ID_134,Intercept]" "r_ID__Hmax[ID_135,Intercept]" "r_ID__Hmax[ID_136,Intercept]"
[46] "r_ID__Hmax[ID_137,Intercept]" "r_ID__Hmax[ID_138,Intercept]" "r_ID__Hmax[ID_139,Intercept]"
[49] "r_ID__Hmax[ID_14,Intercept]" "r_ID__Hmax[ID_142,Intercept]" "r_ID__Hmax[ID_143,Intercept]"
[52] "r_ID__Hmax[ID_144,Intercept]" "r_ID__Hmax[ID_147,Intercept]" "r_ID__Hmax[ID_148,Intercept]"
[55] "r_ID__Hmax[ID_15,Intercept]" "r_ID__Hmax[ID_152,Intercept]" "r_ID__Hmax[ID_153,Intercept]"
[58] "r_ID__Hmax[ID_154,Intercept]" "r_ID__Hmax[ID_155,Intercept]" "r_ID__Hmax[ID_156,Intercept]"
[61] "r_ID__Hmax[ID_157,Intercept]" "r_ID__Hmax[ID_158,Intercept]" "r_ID__Hmax[ID_16,Intercept]"
[64] "r_ID__Hmax[ID_160,Intercept]" "r_ID__Hmax[ID_161,Intercept]" "r_ID__Hmax[ID_163,Intercept]"
[67] "r_ID__Hmax[ID_165,Intercept]" "r_ID__Hmax[ID_166,Intercept]" "r_ID__Hmax[ID_169,Intercept]"
[70] "r_ID__Hmax[ID_17,Intercept]" "r_ID__Hmax[ID_170,Intercept]" "r_ID__Hmax[ID_171,Intercept]"
[73] "r_ID__Hmax[ID_172,Intercept]" "r_ID__Hmax[ID_174,Intercept]" "r_ID__Hmax[ID_175,Intercept]"
[76] "r_ID__Hmax[ID_176,Intercept]" "r_ID__Hmax[ID_177,Intercept]" "r_ID__Hmax[ID_18,Intercept]"
[79] "r_ID__Hmax[ID_180,Intercept]" "r_ID__Hmax[ID_181,Intercept]" "r_ID__Hmax[ID_186,Intercept]"
[82] "r_ID__Hmax[ID_189,Intercept]" "r_ID__Hmax[ID_19,Intercept]" "r_ID__Hmax[ID_191,Intercept]"
[85] "r_ID__Hmax[ID_192,Intercept]" "r_ID__Hmax[ID_197,Intercept]" "r_ID__Hmax[ID_198,Intercept]"
[88] "r_ID__Hmax[ID_199,Intercept]" "r_ID__Hmax[ID_2,Intercept]" "r_ID__Hmax[ID_20,Intercept]"
[91] "r_ID__Hmax[ID_201,Intercept]" "r_ID__Hmax[ID_202,Intercept]" "r_ID__Hmax[ID_207,Intercept]"
[94] "r_ID__Hmax[ID_208,Intercept]" "r_ID__Hmax[ID_21,Intercept]" "r_ID__Hmax[ID_212,Intercept]"
[97] "r_ID__Hmax[ID_214,Intercept]" "r_ID__Hmax[ID_217,Intercept]" "r_ID__Hmax[ID_223,Intercept]"
[100] "r_ID__Hmax[ID_224,Intercept]" "r_ID__Hmax[ID_227,Intercept]" "r_ID__Hmax[ID_23,Intercept]"
[103] "r_ID__Hmax[ID_231,Intercept]" "r_ID__Hmax[ID_234,Intercept]" "r_ID__Hmax[ID_24,Intercept]"
[106] "r_ID__Hmax[ID_25,Intercept]" "r_ID__Hmax[ID_26,Intercept]" "r_ID__Hmax[ID_27,Intercept]"
[109] "r_ID__Hmax[ID_28,Intercept]" "r_ID__Hmax[ID_29,Intercept]" "r_ID__Hmax[ID_3,Intercept]"
[112] "r_ID__Hmax[ID_30,Intercept]" "r_ID__Hmax[ID_32,Intercept]" "r_ID__Hmax[ID_34,Intercept]"
[115] "r_ID__Hmax[ID_35,Intercept]" "r_ID__Hmax[ID_36,Intercept]" "r_ID__Hmax[ID_37,Intercept]"
[118] "r_ID__Hmax[ID_38,Intercept]" "r_ID__Hmax[ID_39,Intercept]" "r_ID__Hmax[ID_4,Intercept]"
[121] "r_ID__Hmax[ID_40,Intercept]" "r_ID__Hmax[ID_41,Intercept]" "r_ID__Hmax[ID_42,Intercept]"
[124] "r_ID__Hmax[ID_45,Intercept]" "r_ID__Hmax[ID_47,Intercept]" "r_ID__Hmax[ID_48,Intercept]"
[127] "r_ID__Hmax[ID_49,Intercept]" "r_ID__Hmax[ID_50,Intercept]" "r_ID__Hmax[ID_51,Intercept]"
[130] "r_ID__Hmax[ID_52,Intercept]" "r_ID__Hmax[ID_53,Intercept]" "r_ID__Hmax[ID_55,Intercept]"
[133] "r_ID__Hmax[ID_56,Intercept]" "r_ID__Hmax[ID_57,Intercept]" "r_ID__Hmax[ID_58,Intercept]"
[136] "r_ID__Hmax[ID_59,Intercept]" "r_ID__Hmax[ID_6,Intercept]" "r_ID__Hmax[ID_60,Intercept]"
[139] "r_ID__Hmax[ID_61,Intercept]" "r_ID__Hmax[ID_62,Intercept]" "r_ID__Hmax[ID_63,Intercept]"
[142] "r_ID__Hmax[ID_64,Intercept]" "r_ID__Hmax[ID_65,Intercept]" "r_ID__Hmax[ID_67,Intercept]"
[145] "r_ID__Hmax[ID_68,Intercept]" "r_ID__Hmax[ID_69,Intercept]" "r_ID__Hmax[ID_7,Intercept]"
[148] "r_ID__Hmax[ID_72,Intercept]" "r_ID__Hmax[ID_73,Intercept]" "r_ID__Hmax[ID_74,Intercept]"
[151] "r_ID__Hmax[ID_75,Intercept]" "r_ID__Hmax[ID_76,Intercept]" "r_ID__Hmax[ID_77,Intercept]"
[154] "r_ID__Hmax[ID_78,Intercept]" "r_ID__Hmax[ID_79,Intercept]" "r_ID__Hmax[ID_8,Intercept]"
[157] "r_ID__Hmax[ID_80,Intercept]" "r_ID__Hmax[ID_81,Intercept]" "r_ID__Hmax[ID_82,Intercept]"
[160] "r_ID__Hmax[ID_84,Intercept]" "r_ID__Hmax[ID_85,Intercept]" "r_ID__Hmax[ID_86,Intercept]"
[163] "r_ID__Hmax[ID_87,Intercept]" "r_ID__Hmax[ID_88,Intercept]" "r_ID__Hmax[ID_89,Intercept]"
[166] "r_ID__Hmax[ID_9,Intercept]" "r_ID__Hmax[ID_90,Intercept]" "r_ID__Hmax[ID_91,Intercept]"
[169] "r_ID__Hmax[ID_92,Intercept]" "r_ID__Hmax[ID_94,Intercept]" "r_ID__Hmax[ID_95,Intercept]"
[172] "r_ID__Hmax[ID_96,Intercept]" "r_ID__Hmax[ID_98,Intercept]" "r_ID__Hmax[ID_99,Intercept]"
[175] "r_ID__k[ID_1,Intercept]" "r_ID__k[ID_10,Intercept]" "r_ID__k[ID_100,Intercept]"
[178] "r_ID__k[ID_101,Intercept]" "r_ID__k[ID_102,Intercept]" "r_ID__k[ID_103,Intercept]"
[181] "r_ID__k[ID_104,Intercept]" "r_ID__k[ID_105,Intercept]" "r_ID__k[ID_106,Intercept]"
[184] "r_ID__k[ID_107,Intercept]" "r_ID__k[ID_109,Intercept]" "r_ID__k[ID_111,Intercept]"
[187] "r_ID__k[ID_112,Intercept]" "r_ID__k[ID_113,Intercept]" "r_ID__k[ID_114,Intercept]"
[190] "r_ID__k[ID_116,Intercept]" "r_ID__k[ID_117,Intercept]" "r_ID__k[ID_118,Intercept]"
[193] "r_ID__k[ID_119,Intercept]" "r_ID__k[ID_12,Intercept]" "r_ID__k[ID_120,Intercept]"
[196] "r_ID__k[ID_121,Intercept]" "r_ID__k[ID_122,Intercept]" "r_ID__k[ID_123,Intercept]"
[199] "r_ID__k[ID_124,Intercept]" "r_ID__k[ID_125,Intercept]" "r_ID__k[ID_126,Intercept]"
[202] "r_ID__k[ID_127,Intercept]" "r_ID__k[ID_128,Intercept]" "r_ID__k[ID_129,Intercept]"
[205] "r_ID__k[ID_13,Intercept]" "r_ID__k[ID_130,Intercept]" "r_ID__k[ID_132,Intercept]"
[208] "r_ID__k[ID_133,Intercept]" "r_ID__k[ID_134,Intercept]" "r_ID__k[ID_135,Intercept]"
[211] "r_ID__k[ID_136,Intercept]" "r_ID__k[ID_137,Intercept]" "r_ID__k[ID_138,Intercept]"
[214] "r_ID__k[ID_139,Intercept]" "r_ID__k[ID_14,Intercept]" "r_ID__k[ID_142,Intercept]"
[217] "r_ID__k[ID_143,Intercept]" "r_ID__k[ID_144,Intercept]" "r_ID__k[ID_147,Intercept]"
[220] "r_ID__k[ID_148,Intercept]" "r_ID__k[ID_15,Intercept]" "r_ID__k[ID_152,Intercept]"
[223] "r_ID__k[ID_153,Intercept]" "r_ID__k[ID_154,Intercept]" "r_ID__k[ID_155,Intercept]"
[226] "r_ID__k[ID_156,Intercept]" "r_ID__k[ID_157,Intercept]" "r_ID__k[ID_158,Intercept]"
[229] "r_ID__k[ID_16,Intercept]" "r_ID__k[ID_160,Intercept]" "r_ID__k[ID_161,Intercept]"
[232] "r_ID__k[ID_163,Intercept]" "r_ID__k[ID_165,Intercept]" "r_ID__k[ID_166,Intercept]"
[235] "r_ID__k[ID_169,Intercept]" "r_ID__k[ID_17,Intercept]" "r_ID__k[ID_170,Intercept]"
[238] "r_ID__k[ID_171,Intercept]" "r_ID__k[ID_172,Intercept]" "r_ID__k[ID_174,Intercept]"
[241] "r_ID__k[ID_175,Intercept]" "r_ID__k[ID_176,Intercept]" "r_ID__k[ID_177,Intercept]"
[244] "r_ID__k[ID_18,Intercept]" "r_ID__k[ID_180,Intercept]" "r_ID__k[ID_181,Intercept]"
[247] "r_ID__k[ID_186,Intercept]" "r_ID__k[ID_189,Intercept]" "r_ID__k[ID_19,Intercept]"
[250] "r_ID__k[ID_191,Intercept]" "r_ID__k[ID_192,Intercept]" "r_ID__k[ID_197,Intercept]"
[253] "r_ID__k[ID_198,Intercept]" "r_ID__k[ID_199,Intercept]" "r_ID__k[ID_2,Intercept]"
[256] "r_ID__k[ID_20,Intercept]" "r_ID__k[ID_201,Intercept]" "r_ID__k[ID_202,Intercept]"
[259] "r_ID__k[ID_207,Intercept]" "r_ID__k[ID_208,Intercept]" "r_ID__k[ID_21,Intercept]"
[262] "r_ID__k[ID_212,Intercept]" "r_ID__k[ID_214,Intercept]" "r_ID__k[ID_217,Intercept]"
[265] "r_ID__k[ID_223,Intercept]" "r_ID__k[ID_224,Intercept]" "r_ID__k[ID_227,Intercept]"
[268] "r_ID__k[ID_23,Intercept]" "r_ID__k[ID_231,Intercept]" "r_ID__k[ID_234,Intercept]"
[271] "r_ID__k[ID_24,Intercept]" "r_ID__k[ID_25,Intercept]" "r_ID__k[ID_26,Intercept]"
[274] "r_ID__k[ID_27,Intercept]" "r_ID__k[ID_28,Intercept]" "r_ID__k[ID_29,Intercept]"
[277] "r_ID__k[ID_3,Intercept]" "r_ID__k[ID_30,Intercept]" "r_ID__k[ID_32,Intercept]"
[280] "r_ID__k[ID_34,Intercept]" "r_ID__k[ID_35,Intercept]" "r_ID__k[ID_36,Intercept]"
[283] "r_ID__k[ID_37,Intercept]" "r_ID__k[ID_38,Intercept]" "r_ID__k[ID_39,Intercept]"
[286] "r_ID__k[ID_4,Intercept]" "r_ID__k[ID_40,Intercept]" "r_ID__k[ID_41,Intercept]"
[289] "r_ID__k[ID_42,Intercept]" "r_ID__k[ID_45,Intercept]" "r_ID__k[ID_47,Intercept]"
[292] "r_ID__k[ID_48,Intercept]" "r_ID__k[ID_49,Intercept]" "r_ID__k[ID_50,Intercept]"
[295] "r_ID__k[ID_51,Intercept]" "r_ID__k[ID_52,Intercept]" "r_ID__k[ID_53,Intercept]"
[298] "r_ID__k[ID_55,Intercept]" "r_ID__k[ID_56,Intercept]" "r_ID__k[ID_57,Intercept]"
[301] "r_ID__k[ID_58,Intercept]" "r_ID__k[ID_59,Intercept]" "r_ID__k[ID_6,Intercept]"
[304] "r_ID__k[ID_60,Intercept]" "r_ID__k[ID_61,Intercept]" "r_ID__k[ID_62,Intercept]"
[307] "r_ID__k[ID_63,Intercept]" "r_ID__k[ID_64,Intercept]" "r_ID__k[ID_65,Intercept]"
[310] "r_ID__k[ID_67,Intercept]" "r_ID__k[ID_68,Intercept]" "r_ID__k[ID_69,Intercept]"
[313] "r_ID__k[ID_7,Intercept]" "r_ID__k[ID_72,Intercept]" "r_ID__k[ID_73,Intercept]"
[316] "r_ID__k[ID_74,Intercept]" "r_ID__k[ID_75,Intercept]" "r_ID__k[ID_76,Intercept]"
[319] "r_ID__k[ID_77,Intercept]" "r_ID__k[ID_78,Intercept]" "r_ID__k[ID_79,Intercept]"
[322] "r_ID__k[ID_8,Intercept]" "r_ID__k[ID_80,Intercept]" "r_ID__k[ID_81,Intercept]"
[325] "r_ID__k[ID_82,Intercept]" "r_ID__k[ID_84,Intercept]" "r_ID__k[ID_85,Intercept]"
[328] "r_ID__k[ID_86,Intercept]" "r_ID__k[ID_87,Intercept]" "r_ID__k[ID_88,Intercept]"
[331] "r_ID__k[ID_89,Intercept]" "r_ID__k[ID_9,Intercept]" "r_ID__k[ID_90,Intercept]"
[334] "r_ID__k[ID_91,Intercept]" "r_ID__k[ID_92,Intercept]" "r_ID__k[ID_94,Intercept]"
[337] "r_ID__k[ID_95,Intercept]" "r_ID__k[ID_96,Intercept]" "r_ID__k[ID_98,Intercept]"
[340] "r_ID__k[ID_99,Intercept]" "r_ID__delta[ID_1,Intercept]" "r_ID__delta[ID_10,Intercept]"
[343] "r_ID__delta[ID_100,Intercept]" "r_ID__delta[ID_101,Intercept]" "r_ID__delta[ID_102,Intercept]"
[346] "r_ID__delta[ID_103,Intercept]" "r_ID__delta[ID_104,Intercept]" "r_ID__delta[ID_105,Intercept]"
[349] "r_ID__delta[ID_106,Intercept]" "r_ID__delta[ID_107,Intercept]" "r_ID__delta[ID_109,Intercept]"
[352] "r_ID__delta[ID_111,Intercept]" "r_ID__delta[ID_112,Intercept]" "r_ID__delta[ID_113,Intercept]"
[355] "r_ID__delta[ID_114,Intercept]" "r_ID__delta[ID_116,Intercept]" "r_ID__delta[ID_117,Intercept]"
[358] "r_ID__delta[ID_118,Intercept]" "r_ID__delta[ID_119,Intercept]" "r_ID__delta[ID_12,Intercept]"
[361] "r_ID__delta[ID_120,Intercept]" "r_ID__delta[ID_121,Intercept]" "r_ID__delta[ID_122,Intercept]"
[364] "r_ID__delta[ID_123,Intercept]" "r_ID__delta[ID_124,Intercept]" "r_ID__delta[ID_125,Intercept]"
[367] "r_ID__delta[ID_126,Intercept]" "r_ID__delta[ID_127,Intercept]" "r_ID__delta[ID_128,Intercept]"
[370] "r_ID__delta[ID_129,Intercept]" "r_ID__delta[ID_13,Intercept]" "r_ID__delta[ID_130,Intercept]"
[373] "r_ID__delta[ID_132,Intercept]" "r_ID__delta[ID_133,Intercept]" "r_ID__delta[ID_134,Intercept]"
[376] "r_ID__delta[ID_135,Intercept]" "r_ID__delta[ID_136,Intercept]" "r_ID__delta[ID_137,Intercept]"
[379] "r_ID__delta[ID_138,Intercept]" "r_ID__delta[ID_139,Intercept]" "r_ID__delta[ID_14,Intercept]"
[382] "r_ID__delta[ID_142,Intercept]" "r_ID__delta[ID_143,Intercept]" "r_ID__delta[ID_144,Intercept]"
[385] "r_ID__delta[ID_147,Intercept]" "r_ID__delta[ID_148,Intercept]" "r_ID__delta[ID_15,Intercept]"
[388] "r_ID__delta[ID_152,Intercept]" "r_ID__delta[ID_153,Intercept]" "r_ID__delta[ID_154,Intercept]"
[391] "r_ID__delta[ID_155,Intercept]" "r_ID__delta[ID_156,Intercept]" "r_ID__delta[ID_157,Intercept]"
[394] "r_ID__delta[ID_158,Intercept]" "r_ID__delta[ID_16,Intercept]" "r_ID__delta[ID_160,Intercept]"
[397] "r_ID__delta[ID_161,Intercept]" "r_ID__delta[ID_163,Intercept]" "r_ID__delta[ID_165,Intercept]"
[400] "r_ID__delta[ID_166,Intercept]" "r_ID__delta[ID_169,Intercept]" "r_ID__delta[ID_17,Intercept]"
[403] "r_ID__delta[ID_170,Intercept]" "r_ID__delta[ID_171,Intercept]" "r_ID__delta[ID_172,Intercept]"
[406] "r_ID__delta[ID_174,Intercept]" "r_ID__delta[ID_175,Intercept]" "r_ID__delta[ID_176,Intercept]"
[409] "r_ID__delta[ID_177,Intercept]" "r_ID__delta[ID_18,Intercept]" "r_ID__delta[ID_180,Intercept]"
[412] "r_ID__delta[ID_181,Intercept]" "r_ID__delta[ID_186,Intercept]" "r_ID__delta[ID_189,Intercept]"
[415] "r_ID__delta[ID_19,Intercept]" "r_ID__delta[ID_191,Intercept]" "r_ID__delta[ID_192,Intercept]"
[418] "r_ID__delta[ID_197,Intercept]" "r_ID__delta[ID_198,Intercept]" "r_ID__delta[ID_199,Intercept]"
[421] "r_ID__delta[ID_2,Intercept]" "r_ID__delta[ID_20,Intercept]" "r_ID__delta[ID_201,Intercept]"
[424] "r_ID__delta[ID_202,Intercept]" "r_ID__delta[ID_207,Intercept]" "r_ID__delta[ID_208,Intercept]"
[427] "r_ID__delta[ID_21,Intercept]" "r_ID__delta[ID_212,Intercept]" "r_ID__delta[ID_214,Intercept]"
[430] "r_ID__delta[ID_217,Intercept]" "r_ID__delta[ID_223,Intercept]" "r_ID__delta[ID_224,Intercept]"
[433] "r_ID__delta[ID_227,Intercept]" "r_ID__delta[ID_23,Intercept]" "r_ID__delta[ID_231,Intercept]"
[436] "r_ID__delta[ID_234,Intercept]" "r_ID__delta[ID_24,Intercept]" "r_ID__delta[ID_25,Intercept]"
[439] "r_ID__delta[ID_26,Intercept]" "r_ID__delta[ID_27,Intercept]" "r_ID__delta[ID_28,Intercept]"
[442] "r_ID__delta[ID_29,Intercept]" "r_ID__delta[ID_3,Intercept]" "r_ID__delta[ID_30,Intercept]"
[445] "r_ID__delta[ID_32,Intercept]" "r_ID__delta[ID_34,Intercept]" "r_ID__delta[ID_35,Intercept]"
[448] "r_ID__delta[ID_36,Intercept]" "r_ID__delta[ID_37,Intercept]" "r_ID__delta[ID_38,Intercept]"
[451] "r_ID__delta[ID_39,Intercept]" "r_ID__delta[ID_4,Intercept]" "r_ID__delta[ID_40,Intercept]"
[454] "r_ID__delta[ID_41,Intercept]" "r_ID__delta[ID_42,Intercept]" "r_ID__delta[ID_45,Intercept]"
[457] "r_ID__delta[ID_47,Intercept]" "r_ID__delta[ID_48,Intercept]" "r_ID__delta[ID_49,Intercept]"
[460] "r_ID__delta[ID_50,Intercept]" "r_ID__delta[ID_51,Intercept]" "r_ID__delta[ID_52,Intercept]"
[463] "r_ID__delta[ID_53,Intercept]" "r_ID__delta[ID_55,Intercept]" "r_ID__delta[ID_56,Intercept]"
[466] "r_ID__delta[ID_57,Intercept]" "r_ID__delta[ID_58,Intercept]" "r_ID__delta[ID_59,Intercept]"
[469] "r_ID__delta[ID_6,Intercept]" "r_ID__delta[ID_60,Intercept]" "r_ID__delta[ID_61,Intercept]"
[472] "r_ID__delta[ID_62,Intercept]" "r_ID__delta[ID_63,Intercept]" "r_ID__delta[ID_64,Intercept]"
[475] "r_ID__delta[ID_65,Intercept]" "r_ID__delta[ID_67,Intercept]" "r_ID__delta[ID_68,Intercept]"
[478] "r_ID__delta[ID_69,Intercept]" "r_ID__delta[ID_7,Intercept]" "r_ID__delta[ID_72,Intercept]"
[481] "r_ID__delta[ID_73,Intercept]" "r_ID__delta[ID_74,Intercept]" "r_ID__delta[ID_75,Intercept]"
[484] "r_ID__delta[ID_76,Intercept]" "r_ID__delta[ID_77,Intercept]" "r_ID__delta[ID_78,Intercept]"
[487] "r_ID__delta[ID_79,Intercept]" "r_ID__delta[ID_8,Intercept]" "r_ID__delta[ID_80,Intercept]"
[490] "r_ID__delta[ID_81,Intercept]" "r_ID__delta[ID_82,Intercept]" "r_ID__delta[ID_84,Intercept]"
[493] "r_ID__delta[ID_85,Intercept]" "r_ID__delta[ID_86,Intercept]" "r_ID__delta[ID_87,Intercept]"
[496] "r_ID__delta[ID_88,Intercept]" "r_ID__delta[ID_89,Intercept]" "r_ID__delta[ID_9,Intercept]"
[499] "r_ID__delta[ID_90,Intercept]" "r_ID__delta[ID_91,Intercept]" "r_ID__delta[ID_92,Intercept]"
[502] "r_ID__delta[ID_94,Intercept]" "r_ID__delta[ID_95,Intercept]" "r_ID__delta[ID_96,Intercept]"
[505] "r_ID__delta[ID_98,Intercept]" "r_ID__delta[ID_99,Intercept]" "lp__"
r_ID__Hmax
datalist0 = list()
for (samples in sample.list) {
ttt0 <- paste0("r_ID__Hmax[",samples,",Intercept]")
dat0 <- extract(width.best.fitted, ttt0, permuted=TRUE)
datalist0[[ttt0]] <- dat0
}
list.data0 = do.call(rbind, datalist0)
Hmax.ID <- data.frame(matrix(unlist(list.data0), nrow = 166, byrow = T))
row.names(Hmax.ID) <- paste0("r_ID__Hmax[",sample.list,",Intercept]")
head(Hmax.ID)
r_ID__k
datalist1 = list()
for (samples in sample.list) {
ttt1 <- paste0("r_ID__k[",samples,",Intercept]")
dat1 <- extract(width.best.fitted, ttt1, permuted=TRUE)
datalist1[[ttt1]] <- dat1
}
list.data1 = do.call(rbind, datalist1)
k.ID <- data.frame(matrix(unlist(list.data1), nrow = 166, byrow = T))
row.names(k.ID) <- paste0("r_ID__k[",sample.list,",Intercept]")
head(k.ID)
r_ID__delta
datalist2 = list()
for (samples in sample.list) {
ttt2 <- paste0("r_ID__delta[",samples,",Intercept]")
dat2 <- extract(width.best.fitted, ttt2, permuted=TRUE)
datalist2[[ttt2]] <- dat2
}
list.data2 = do.call(rbind, datalist2)
delta.ID <- data.frame(matrix(unlist(list.data2), nrow = 166, byrow = T))
row.names(delta.ID) <- paste0("r_ID__delta[",sample.list,",Intercept]")
head(delta.ID)
estimates of parameters
par.list = c("b_Hmax_Intercept", "b_k_Intercept", "b_delta_Intercept", "b_Hmin_Intercept", "sd_ID__Hmax_Intercept", "sd_ID__k_Intercept","sd_ID__delta_Intercept","sigma", "lp__")
datalist3 = list()
for (pars in par.list) {
dat3 <- extract(width.best.fitted, pars, permuted=TRUE)
datalist3[[pars]] <- dat3
}
list.data3 = do.call(rbind, datalist3)
pars.est <- data.frame(matrix(unlist(list.data3), nrow = 9, byrow = T))
row.names(pars.est) <- par.list
pars.est
combine estimates of all parameters
pars.estimate.width.fitted <- rbind(pars.est, Hmax.ID, k.ID, delta.ID)
head(pars.estimate.width.fitted, 10)
write.csv(pars.estimate.width.fitted, file="pars.estimate.width.fitted.csv")